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ABSTRACT 

Outward migration of low-mass planets has recently been shown to be a possibility in non-barotropic disks. We 
examine the consequences of this result in evolutionary models of protoplanetary disks. Planet migration occurs 
toward equilibrium radii with zero torque. These radii themselves migrate inwards because of viscous accretion and 
photoevaporation. We show that as the surface density and temperature fall the planet orbital migration and disk 
depletion timescales eventually become comparable, with the precise timing depending on the mass of the planet, 
when this occurs, the planet decouples from the equilibrium radius. At this time, however, the gas surface density 
is already too low to drive substantial further migration. A higher mass planet, of 10 Mqj, can open a gap during the 
late evolution of the disk, and stops migrating. Low-mass planets, with 1 or 0.1 Mgj , released beyond 1 AU in our 
models avoid migrating into the star. Our results provide support for the reduced migration rates adopted in recent 
planet population S5mthesis models. 

Subject headings: 



1. INTRODUCTION 

The origin of planetary systems remains a major chal- 
lenge to astrophysical theory. Aside from the quandary 
of planet formation, planet survival is also a prob- 
lem. Planet formation occurs when the gas disk is still 
present, and by exchanging angular momentum with 
the gas, planets start to migrate. This occurs as the 
planet excites perturbations in the disk that, in turn, ex- 
ert torques on the planet. The asymmetry of these per- 
turbations on either side of the planet determines the 
strength of the torques, and thus the direction of migra- 
tion. 

The main components of the excited perturbations are 
the one-armed spirals launched at the Lrndblad reso- 
nances and the librating material on horseshoe orbits 
in the planet's corotation region. For low-mass plan- 
ets, with weak wakes, both the shape of the spiral wake 
and the resulting Lindblad torque can be treated by lin- 
ear analysis under some assumptions, particularly of 
local isothermality. The analytical prediction (Goldre- 
ich & Tremaine 1979; Ward 1997; Tanaka et al. 2002), 
confirmed by numerical simulations (Nelson et al. 2000; 
D'Angelo et al. 2002; Bate et al. 2003), is that the inner 
Lindblad resonances lead to positive torques, whereas 
the outer ones lead to negative torques. The outer Lind- 
blad resonances lie closer to the planet, and thus pro- 
duce stronger torques, so the planet migrates inwards. 
This migration mode, referred to as Type I migration, 

occurs on timescales between 10^ and 10^ yr. This 
is a serious problem for planet formation since these 
timescales are much shorter than the lifetimes of disks 
(10^-10'^ yr). Halting or slowing Type I migration is im- 
perative if planets are to survive at all. Indeed, planet 
population synthesis models (e.g., Alibert et al. 2004; 
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Mordasini et al. 2009) have to assume a reduction fac- 
tor for Type I migration of 30-1000 in order to match the 
observed distribution of planetary semimajor axes. 

Paardekooper & Mellema (2006) made a major step 
toward the solution of this problem. They found 
that when the locally isothermal approximation usu- 
ally assumed in the literature was relaxed, the plan- 
ets migrated outwards. This behavior was explained 
by Baruteau & Masset (2008) and Paardekooper & Pa- 
paloizou (2008) as resulting from an entropy-related 
torque exerted by material on horseshoe orbits in the 
corotation region. This mechanism operates in regions 
of the disk that have a negative entropy gradient and 
inefficient radiative cooling, where sustenance of the 
torque requires some viscous and thermal diffusion 
(Paardekooper & Papaloizou 2008). Kley & Crida (2008) 
and Kley et al. (2009) showed that outward migration 
indeed occurs in disks with realistic heating and cool- 
ing. 

This outward migration remains rapid. Planets mi- 
grate toward the outer disk, reach an equilibrium radius 
of zero torque, and stay put thereafter. As emphasized 
by Paardekooper et al. (2010), slow inward migration 
then occurs as disk evolution shifts the equilibrium ra- 
dius inwards. The situation becomes similar to that of 
gap-opening planets (Type II migration), with the planet 
migrating in lockstep with the gas as the gas accretes. 

This scenario poses a subtle problem that we examine 
in this Letter. If the planet is moving together with the 
gas as the disk depletes, at some point in the evolution 
the disk may reach a thermodynamic state such that in- 
ward migration resumes. This brings the problem back 
to square one, because, if the planet is to survive, the 
remaining disk lifetime must be shorter than the Type 
I migration timescale in that evolutionary state. We ex- 
amine this possibility using one-dimensional evolution- 
ary models of protoplanetary disks including heating 
and cooling. We describe the model in the following 
section, present our results in § 3, and give concluding 
remarks in § 4. 
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2. THE MODEL 

2.1. Gas evolution 

We consider non-irradiated disks evolving by viscous 
diffusion and photoevaporation (Lynden-Bell & Pringle 
1974; Lin & Papaloizou 1979), 
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where S is the surface density, r is the radius, and v 
is the effective viscosity. We take the photoevaporation 
rate to be (Veras & Armitage 2004; Mordasini et al. 2009) 

0, for r < Rg, 

My,/[2n{rext-Rg)r], for r > J?^, 

which is valid for external irradiation. 

For temperature evolution, we use a model without 
shock heating (Nakamoto & Nakagawa 1994) 
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where T and T;, are the midplane and background tem- 
peratures, respectively, Q is the Keplerian frequency, 
and (7 is the Stefan-Boltzmann constant. We take the 
effective optical depth at the midplane (Hubeny 1990; 
Kley & Crida 2008): 
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Equation (3) states that the emerging flux is the result 
of an equilibrium between viscous heating, background 
radiation, and radiative cooling. The optical depth is 
T = kS /2, and the opacities k are taken from Bell et 
al. (1997). We assume that although dust growth and 
planet formation lock away refractory material, frag- 
mentation efficiently replenishes small grains, keeping 
the disks opaque during their evolution (Birnstiel et al. 
2009). 

2.2. Planet evolution 
The planet's orbital radius evolves as 
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where nip is the planet's mass and F is the torque 
from the gas. We assume a circular orbit and constant 
planet mass. We only consider low-mass planets, and 
so can ignore the back reaction of the planet onto the 
gas (Alexander et al. 2006). The torques are modeled 
with analytical fits valid for the fully imsaturated case 
(Paardekooper et al. 2010). Using b/h=OA, where h is 
the disk's aspect ratio and b is the gravitational soften- 
ing of the planet's potential (necessary in models with 
less than three dimensions), the torques are 

riso/ro = -0.85-a-0.9^, (6) 
7rad/ro =-0.85- a -1.7/5 + 7.9^/7 (7) 

for the locally isothermal and adiabatic equations of 
state, respectively. The adiabatic index 7 = 1.4, and a. 



/5, and ^ are the negative of the local density, tempera- 
ture, and entropy gradients: 
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The torques are normalized by 



^ = /3-(7-l)a. (8) 
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where c\ is the ratio of planetary to stellar mass and Z'p 
is the surface density at the position of the planet. We 
interpolate between the two torque regimes to get 
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are the ra- 



Here, = frad/^dyn/ where t^^^ and tjyn 
diative and dynamical timescales. If 6* -C 1, radia- 
tive processes can restore the temperature quickly com- 
pared to the horseshoe turnover time, leading to isother- 
mal horseshoe turns. It was shown in Paardekooper 
& Papaloizou (2009) that the timescale to establish the 
full horseshoe drag (fdrag) is a fraction of the libration 
timescale (typically 10 tdyn for a 5 planet). How- 
ever, as long as fj-ad > ^dyn/ the torque is affected by non- 
linear effects, and since the nonlinear torque is so much 
stronger than its linear counterpart (Paardekooper et al. 
2010), the parameter governing linearity (and therefore 
isothermality) is Q rather than frad/^drag- 

Since f^rag depends on c\, using f^rag would lead to 
different migration behavior for different planet masses. 
In practice, the transition between isothermal and adia- 
batic regions occurs on such a small length scale that this 
would not change the qualitative outcome of the mod- 
els. Mass segregation is more likely to come from effects 
of saturation, which we do not consider in this Letter, 
since it is as yet poorly understood for non-barotropic 
disks. We note, however, that for viscously heated, op- 
tically thick disks close to thermal equilibrium, radia- 
tive thermal diffusion will approach the viscous heating 
rate. Therefore, even in the very optically thick inner re- 
gions of the disk, a sizeable fraction of the unsaturated 
torque can be sustained as long as the disk remains vis- 
cous. 

To calculate 0, we consider f j-ad = where £ is the 
internal energy. The cooling E = V ■ F, where F is the 
flux. Using E = c^pT , p = S/2H, and |F 
have 



o-T^ii, we 



^rad 

The dynamical time t^. 
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2.3. Simulation parameters 

We use a one-dimensional linear grid covering 0.1- 
30 AU with 200 points. The surface density is specified 
by the initial mass accretion rate Mq, and the viscosity 
parameter ass (Shakura & Sunyaev 1973), following the 
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analytical fits of Papaloizou & Terquem (1999) for S- 
V relations. We use Mq = 10"'' yr"^ and ass = 10"^. 
The wind is modeled with Mw = 10~^ M0 yr~^ and = 
5AU. For a given surface density. Equation (3) speci- 
fies the temperature, with Tf, = 10 K. Because the optical 
depth depends on temperature, we solve Equation (3) 
with a Newton-Raphson root-finding algorithm (using 
0.01 K precision). We examine planets of mass 0.1, 1, and 
10 Mg). Planet-planet interactions are ignored. Bound- 
ary conditions are taken as outflow. We compute the 
derivatives as in the PENCIL CODE^, with sixth-order 
spatial derivatives and a third-order Rimge-Kutta time 
integrator. 

3. RESULTS 

Figure 1 shows the evolution of our fiducial disk 
model. The changes in density slope initially located at 
~0.4 and 21 AU correspond to opacity transitions (Pa- 
paloizou & Terquem 1999). The constant temperature 
plateau initially at 4-5 AU corresponds to the opacity 
transition at ~130 K where ice grains sublimate (see Fig- 
ure Ic). Initially, the ratio of radiative to dynamical 
timescales 0>1 through most of the disk (6* < 1 only at 
r > 70 AU, beyond the model grid). Thus, the torque in 
Equation (10) is mostly adiabatic, and largely indepen- 
dent of the interpolation procedure. 

Total disk mass and accretion rate decay nearly lin- 
early with time. A disk of 0.08 Mr^, initially accreting at 
10"^ Mq yr"^ gets depleted in 4.8 Myr. At that time, the 
total disk mass is < 1 M®, and we terminate the simu- 
lation. The model gives a power-law decline of the den- 
sity. An inner hole (Clarke et al. 2001) never forms, be- 
cause of the shallow radial slope of unity for the wind 
driven by external photoevaporation. 

As the surface density drops, so does the viscous 
heating rate (Equation 3), and the temperature falls ac- 
cordingly, down to T = Tj,. During the evolution, the 
isothermal plateau shifts inwards from 4 AU to 1 AU by 
4.5 Myr. The optical depth only drops to unity at 10 AU 
after 4 Myr, and at 5 AU after 4.4 Myr, so the planet for- 
mation region is optically thick through most of the disk 
evolution. To check our method, we also ran models 
with Equation (3), and models where the temperature 
was derived using (l+l)-dimensional models, follow- 
ing Papaloizou & Terquem (1999). The differences were 
minor, consisting of a slightly larger isothermal plateau 
(extending to 6 AU instead of 5 AU), and an outer disk 
~10K hotter. 

To understand our results for planet migration in 
evolving disks, we first examine the behavior of planets 
in stationary models. In Figure 2, we show the torque 
as a function of radius. Migration halts at stable equilib- 
rium radii where F = and dV / dr < 0, that is, where the 
torque is positive within (corresponding to outward mi- 
gration) and negative outside (inward migration). Two 
such equilibrium radii occur, at 4AU and 21 AU, cor- 
responding to the inner boundaries of the isothermal 
plateau and the outer isothermal region, where T = T;,. 
(There is also a thin region around an equilibrium ra- 
dius at 0.6 AU that we do not consider here.) A negative 
torque acts on planets migrating from Yp > 21 AU, bring- 

^ See http: / / www.nordita.org/ software /pencil-code 



ing them to that radius, while a positive torque acts on 
planets with 4 < fp / AU < 21, also bringing them to that 
radius. Planets with fp < 4 AU migrate to the inner equi- 
librium radius. 

We released planets of varying mass at several radii in 
the disk, and see precisely this behavior (Figure 3a). The 
migration time to the equilibrium radius is inversely 
proportional to the planet mass, because the torque de- 
pends quadratically on mass (Equation 9), so Tp <x l/nip 
(Equation 5). Planets of 10 reach the equilibrium 
radii in < 0.1 Myr. In the isothermal case, planets mi- 
grate inwards at comparable speeds. 

In Figure 3b, we show trajectories of planet migra- 
tion in an evolving disk. Now the radii of zero torque 
shift inwards as the disk evolves. The planets migrate 
to these radii on the Type I timescale of 10^ yr, and then 
couple to disk evolution. Subsequently, they migrate in- 
wards on the slow accretion timescale of 10^ yr, compa- 
rable to Type II migration, as predicted by Paardekooper 
et al. (2010). 

Figure 3b has some intriguing features that give 
pause. There is a clear effect of mass on the later evo- 
lution of the planets. In the stationary case, the only 
effect of mass is to determine the speed of migration to 
the equilibrium radii. In the evolving case, the mass also 
helps determine the final location of the planet. To show 
this, we separately plot the tracks of planets of different 
masses in Figures 3c-e, as well as the trajectories of the 
equilibrium radii. The planets follow these equilibrium 
radii until late times, when they decouple. 

We can understand the process of decoupling by com- 
paring the viscous accretion timescale to the migration 
timescale. As the disk thins, the surface density reaches 
a value so low that the gas cannot transfer sufficient an- 
gular momentum to the planet for its orbital radius to 
evolve as fast the equilibrium radius. Another way of 
understanding this is that a perturbation of the planet 
away from the equilibrium radius will only be corrected 
if the torque is sufficiently strong to return the planet 
before the radius moves a substantial distance. As the 
torque is a function of planet mass, the time and radius 
when decoupling occurs are also functions of the planet 
mass. 

We can estimate this time and radius by comparing 
the migration timescale {tmig=^p/\^p\) to the disk accre- 
tion timescale {ty=S As long as tmig < tv, the 
planet can keep up with the evolution of the disk. When 
this is no longer true, decoupling occurs, releasing the 
planet. We plot the quantity tmig/^i' in Figure 4, for the 
planet of 1 Mfg at different times. The structure in the 
plot comes from the radial derivatives of Sv that de- 
fine the mass accretion flow, and from the torques F that 
define the migration rate. The spikes in the figure cor- 
respond to the equilibrium locations, where F = (and 
therefore formally tmig/ iu = The decoupling in Fig- 
ure 3d starts to occur at w 4.0 Myr. In Figure 4, we see 
that it roughly corresponds to the time when the vicinity 
of the equilibrium radius crosses the line of ^mig/ ty = i- 

Evolution after decoupling proceeds as follows. The 
inner equilibrium radius corresponds to the inner edge 
of the isothermal plateau. Once the planet decouples, it 
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Fig. 1. — Evolution of (a) surface density, (b) midplane temperature, (c) optical depth, and (d) 0, the ratio of radiative to dynamical timescales. 
The disk is drained after 4.8 Myr, with mass dropping below Ma. The temperature shows a plateau at »il30 K, produced by the opacity transition 
when ice grains sublimate that moves inward as the disk cools. Through most of the evolution, > 1, so torques are adiabatic. 




r{AU) 

Fig. 2. — Radial distribution of torque over time (Equation 10). Pos- 
itive torques drive outward migration. The equilibrium radii of zero 
torque shift inwards as the disk evolves. The isothermal plateau cor- 
responds to the well of negative torques initially at 4—5 AU. Discon- 
tinuities correspond to opacity jumps. Torques are scaled by IJ/T to 
aid visualization. 

is released inside the plateau, itself a region of inward 
migration since the temperature gradient vanishes. As 
disk evolution continues, the decoupled planet soon 
finds itself at the outer edge of the isothermal plateau, 
and starts migrating outwards. For planets of 1 Mfy, de- 
coupling occurs at r^=l AU at 4.5 Myr (Figure 3b). The 
planet then rapidly descends the temperature gradient 
until it reaches the outer equilibrium radius. However, 
this radius too moves inward faster than the planet can 
migrate, so the planet enters the isothermal outer disk 



(where T ~ Tj,), another region of slow inward migra- 
tion. At this stage, a planet of 1 Mg^ still has time to mi- 
grate from Vp = 2.0 AU to Vp = 1.5 AU before the density 
drops too low to drive further migration (Figure 3d). 

The planet of 10 is strongly coupled and would 
follow the outer equilibrium radius until it hit the inner 
boundary of our model at 0.1 AU. The halt at 2 AU seen 
in Figure 3e represents an artificial termination of the 
simulation. At that point, two criteria had been fulfilled. 
First, the scale height had become smaller than the Hill 
radius of the planet, so gap formation should have oc- 
curred. This is not sufficient to terminate the simula- 
tion, because the inward motion of the equilibrium ra- 
dius itself occurs at the timescale of Type 11 migration. 
However, the second criterion was that the mass param- 
eter determining Type 11 migration, Sr^ (Mordasini et 
al. 2009), had become smaller than the planet's mass. At 
this stage. Type 11 migration becomes planet dominated 
and we consider that it comes to a halt. The smaller 
planets never carve gaps. 

We investigate migration in disks with different val- 
ues of Mw and Mq, yet constrained by a lifetime of 
1-10 Myr, without finding qualitative differences. A 
change in migration behavior is only seen for hotter 
disks. For ass=0.1, the disks show 6* w 1 throughout. 
Migration then shows a mix of isothermal and adiabatic 
behavior, being mostly inwards yet with equilibrium 
points present. The torque only becomes isothermally 
dominated for ass —1/ which yields an unrealistically 
high accretion rate. 

4. CONCLUSIONS 
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Fig. 3. — (a) Orbital migration in stationary disks. Planets migrate at mass-dependent rates toward equilibrium radii of zero torque (see Figure 2). 
(b) Orbital migration in evolving disks. Equilibrium radii migrate inwards on the accretion timescale, but the planets eventually migrate too slowly 
to remain in equilibrium. In the lower panels (c)-(e) we show the same as in panel (b) but separated by planet mass. The tracks of the equilibrium 
radii (dotted lines) are shown for comparison. 
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Fig. 4. — Comparison of the timescales of planet migration (tmig) and 
disk accretion (fi,) for a planet of 1 M®. 



In this Letter, we examine the trajectories of plan- 
ets undergoing Type 1 migration in evolving, radiative 
disks with initially nearly adiabatic midplanes, where 
outward migration can occur. Planets migrate toward 
equilibrium radii, where the torque acting on them van- 
ishes. These radii correspond to opacity jumps and 
to the transition, in the outer disk, to an isothermal 
state. Because of viscous accretion and photoevapora- 
tion, these equilibrium radii themselves move inwards 
on disk accretion timescales. 

As long as the torques are strong enough to keep the 



planets coupled to disk evolution, the planets migrate 
in lockstep with the gas at the accretion timescale. How- 
ever, as the disk surface density drops, the timescales of 
orbital migration and disk accretion eventually become 
comparable. At this stage, if the planet is perturbed 
from an inner equilibrium radius, the equilibrium ra- 
dius moves inward faster than the torques can return 
the planet, so it decouples. 

If the continuing migration is outwards, the planet en- 
counters another equilibrium radius, and the same pro- 
cess of locking and decoupling occurs. The outermost 
equilibrium radius lies where the disk reaches the back- 
ground temperature, and from there the planet can only 
migrate inwards. However, by the time that the planet 
decouples even from that equilibrium radius, the disk 
is already so severely depleted that the ensuing inward 
migration is feeble, and soon comes to a halt, as the re- 
maining disk mass cannot exert a substantial torque. In 
no case did a planet released beyond 1 AU migrate all 
the way to the star. 

We stress that we only use a single model for the opac- 
ities, which may change as planet formation progresses. 
The dust size distribution must depend on the balance 
between coagulation and fragmentation, which remains 
poorly understood. We also neglect stellar irradiation, 
which will become important in the late stages of the 
disk evolution. Irradiation will maintain high temper- 
atures in the inner disk (~100K; Chambers 2009), pre- 
venting gap formation for 10 planets. On the other 
hand, irradiation should lead to a hole-forming photo- 
evaporating wind (Clarke et al. 2001; Alexander et al. 
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2006), which quickly depletes the disk, thus possibly 
bringing migration to an even earlier halt. Future work 
should self -consistently address these issues. 

As a consequence of the independence of equilibrium 
radii on planet mass, all planets migrate to these equi- 
librium locations. Ensembles of planets reaching them 
may become violently unstable due to N-body inter- 
actions. Nevertheless, even if scattered away, migra- 
tion will invariably drive the planets back toward these 
radii. The final outcome may well be collisions driving 
further planet growth, aiding rapid giant planet forma- 
tion or forming planets in 1:1 resonance. 

If this is the case, however, it raises the question of 
why the solar system has a set of neatly spaced plan- 
ets as opposed to only two, as the two equilibrium radii 
of the model might naively suggest. One possible solu- 
tion is that r (and thus any equilibrium radius) shows 
a dependence on the planet-to-star mass ratio q at the 
verge of gap opening when q ^ (Masset et al. 2006). 
Another is that we only consider the fully unsaturated 
torque, whereas saturation depends on the width of the 
horseshoe region and therefore on the planet's mass. 
However, the level of saturation in radiative disks is not 
fully understood at present, and we cannot easily add 
it to our study. Future models should include effects of 
saturation to study possible mass segregation. Finally, 
in view of the long migration timescales for M < 0.1 Mg^, 



such planets may just not have the time to migrate back 
to the equilibrium location before the disk vanishes if 
scattered far enough. This scattered population of small 
planets could provide the initial conditions for the ter- 
restrial planets of our own solar system. 

Our results provide qualitative and quantitative jus- 
tification for the reduction of Type I migration rates as- 
sumed in planetary population synthesis models (e.g., 
Alibert et al. 2004; Ida & Lin 2008; Mordasini et al. 
2009). Instead of migrating on the fast, mass-dependent, 
timescale ^mig/ we find that planets spend their first 
Myr near equilibrium radii that change only on the 
slow accretion timescale ty. We show in Figure 4 that 
fmig/ ty^O.l for a 1 planet during most of the evolu- 
tion of the disk. Examining the same figure for different 
masses shows a linear dependence on mass, tmig / tp ^ 
0.1 (M/ M®), consistent with the population synthesis 
assumptions. 
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